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ABSTRACT 

We present a detailed non-spherical modeling of dark matter halos on the 
basis of a combined analysis of the high- resolution halo simulations (12 halos 
with N ~ 10 6 particles within their virial radius) and the large cosmological sim- 
ulations (5 realizations with N = 512 3 particles in a 100/i _1 Mpc boxsize). The 
density profiles of those simulated halos are well approximated by a sequence of 
the concentric triaxial distribution with their axis directions being fairly aligned. 
We characterize the triaxial model quantitatively by generalizing the universal 
density profile which has previously been discussed only in the framework of the 
spherical model. We obtain a series of practically useful fitting formulae in ap- 
plying the triaxial model; the mass and redshift dependence of the axis ratio, the 
mean of the concentration parameter, and the probability distribution functions 
of the the axis ratio and the concentration parameter. These accurate fitting 
formulae form a complete description of the triaxial density profiles of halos in 
Cold Dark Matter models. Our current description of the dark halos will be 
particularly useful in predicting a variety of nonsphericity effects, to a reason- 
ably reliable degree, including the weak and strong lens statistics, the orbital 
evolution of galactic satellites and triaxiality of galactic halos, and the non-linear 
clustering of dark matter. In addition, this provides a useful framework for the 
non-spherical modeling of the intra-cluster gas, which is crucial in discussing the 
gas and temperature profiles of X-ray clusters and the Hubble constant estimated 
via the Sunyaev - Zel'dovich effect. 
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Subject headings: galaxies: clusters: general - cosmology: miscellaneous - meth- 
ods: numerical 



1. Introduction 

The density profiles of dark matter halos have attracted a lot of attention recently since 
Navarro, Frenk, & White (1996, 1997, NFW hereafter) discovered the unexpected scaling 
behavior in their simulated halos. Subsequent independent higher-resolution simulations 
(e.g., Fukushige & Makino 1997, 2001; Moore et al. 1998; Jing 2000; Jing & Suto 2000) 
confirmed the validity of the NFW modeling, in particular the presence of the central cusp, 
although the inner slope of the cusp seems somewhat steeper than they originally claimed. 
Those previous models, however, have been based on the spherical average of the density 
profiles. Actually it is also surprising that the fairly accurate scaling relation applies after 
the spherical average despite the fact that the departure from the spherical symmetry is 
quite visible in almost all simulated halos (e.g., Fig.l of Jing & Suto 2000). 

A more realistic modeling of dark matter halos beyond the spherical approximation is 
important in understanding various observed properties of galaxy clusters and non-linear 
clustering (especially the high-order clustering statistics) of dark matter in general. In par- 
ticular, the non-sphericity of dark halos is supposed to play a central role in the X-ray 
morphologies of clusters (Jing, et al. 1995; Buote & Xu 1997), in the cosmological param- 
eter determination via the Sunyaev-Zel'dovich effect (Birkinshaw, Hughes & Arnaud 1991; 
Inagaki, Suginohara & Suto 1995; Yoshikawa, Itoh & Suto 1998) and in the prediction of the 
cluster weak lensing and the gravitational arc statistics(Bartelmann et al. 1998; Meneghetti 
et al. 2000, 2001; Molikawa & Hattori 2001; Oguri, Taruya, & Suto 2001; Keeton & Madau 
2001). Nevertheless useful analytical modeling of the non-sphericity is almost impossible, 
and numerical simulations are the only practical means to provide statistical information. 

While the non-sphericity of the dark matter halos is a poorly studied topic, some seminal 
studies do exist which attempt to detect and characterize the non-spherical signature (e.g., 
Barnes & Efstathiou 1987; Warren et al. 1992; Dubinski 1994; Jing, et al. 1995; Thomas et 
al. 1998; Yoshida et al. 2000; Meneghetti et al. 2001; Bullock 2001). Nevertheless there is no 
systematic and statistical study to model and characterize the density profiles of simulated 
halos. This is exactly what we will present in the rest of the paper. In particular, much higher 
mass and spatial resolutions of our current N-body simulations enable us to characterize the 
statistics of the halo non-sphericity with an unprecedented precision. 

This paper is organized as follows; two different sets of N-body simulations that we 
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extensively analyze here are described in §2. In §3, we discuss how to define the iso-density 
surfaces of dark mater halos from simulation data and then argue that they are well approx- 
imated by a sequence of the concentric triaxial model. Section 4 characterizes the statistical 
distribution of the triaxial model parameters. Finally §5 is devoted to summary and discus- 
sion. 



2. Simulations for dark matter halos 

We use two different simulations for the current purpose. The first is our new set of 
cosmological N-body simulations with iV = 512 3 particles in a 100/i _1 Mpc box, and the other 
is a set of high-resolution halo simulation runs. We describe the two simulations in the next 
subsections in order. 



2.1. Cosmological simulations 

The first set of simulations is our new runs with N = 512 3 particles in a 100/i _1 Mpc box. 
These runs have been carried out in 2001 with our Particle-Particle-Particle-Mesh (P 3 M) code 
on the vector-parallel machine VPP5000 at the National Astronomical Observatory of Japan. 
The code adopts the standard P 3 M algorithm (Hockney & Eastwood 1981; Efstathiou et 
al. 1985), is vectorized (Jing & Suto 1998), and has been recently parallelized. A mesh of 
1200 3 grid points is used for the Particle-Mesh (PM) force computation with the optimized 
Green function (Hockney & Eastwood 1981). The short-range force is compensated for the 
PM force calculation at the separation less than e = 2.7 H, where H is the mesh cell size 
(Efstathiou et al. 1985). The linked-list technique has been used for computing the short- 
range Particle-Particle (PP) interaction with 448 3 linked-list cells. The computer has a total 
of 64 processors, and we use iVcpu = 8 to 32 processors, upon their availability, to run our 
code. The most important advantage of the machine for our work is that each processor 
has a big memory of 16GB, sufficient for storing all the information of the code. The PM 
computation can be easily parallelized, and it is crucial to parallelize the PP computation 
which dominates the CPU computation time for a strongly clustered simulation like our 
present case. We sliced the simulation box in one direction (e.g. z-axis) with the thickness 
chosen to be the cell size of the linked-list cell. Those 448 slices in total are sorted in the 
descending order according to the number of particles they contain. We distribute the PP 
force computation among the different processors in a simple way; n-th processor {n — 1 to 
iVcpu) is assigned the force computation for those slices with indices of jNcpv+n where j runs 
from to j max = 448/iVcpu — 1- The PP interaction of the particles in the same slice and in 
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the adjacent lower slice is considered, so the interaction for each pair of particles is computed 
only once. With this computation partition, we find that the load-balance problem, which 
becomes progressively serious for P 3 M simulations in the later strongly clustered regime, can 
be overcome to a satisfactory degree; even at the final stage of our simulation runs, the CPU 
time for the PP part is nearly inversely proportional to the number of the processors used. 
This implies that the code has achieved a good parallelization efficiency 

We consider two representative cold dark matter (CDM) models; a low-density flat 
cosmological model (LCDM) with Q = 0.3 and Ao = 0.7, and the Einstein-de Sitter model 
with Qq = 1 (SCDM). The primordial density fluctuation is assumed to obey the Gaussian 
statistics, and the power spectrum is given by the Harrison-Zel'dovich type. The linear 
transfer function for the dark matter power spectrum is taken from Bardeen et al. (1986). 
The shape and the normalization of the linear power spectrum are specified by the shape 
parameter, T = £l h, and a 8 respectively, where h is the Hubble constant in lOOkms^Mpc -1 
and cr 8 is the rms linear density fluctuation within the sphere of the radius 8 /i _1 Mpc. Table 1 
summarizes the physical and simulation parameters used for these simulations. We adopted 
eg = 0.9 for LCDM and 0.55 for SCDM, both of which are slightly smaller than those in our 
previous simulations (Jing & Suto 1998), but seem more consistent with recent observations 
(e.g., Seljak 2002; Lahav et al. 2002). With the adopted values for those physical parameters, 
the LCDM model satisfies almost all current observations while the SCDM model is known 
to have many difficulties. Therefore we mainly analyze the LCDM model for our purpose, 
and sometimes use the SCDM simulation just for comparison. 

The boxsize of our cosmological simulations is 100 /i _1 Mpc, so the particle mass is m p = 
6.2 x 10 8 /i- 1 M Q and 2.1 x lO 9 /i _1 M , respectively, for the LCDM and SCDM simulations 
(Table 1). The force resolution is rj = 20 h^kpc for the linear density softening form 
(Efstathiou et al. 1985; this roughly corresponds to rj/3 for the Plummer-type softening 
length). The simulations are evolved by 1200 time steps from the initial redshift Z; L = 72. 
Two realizations are computed for each model. One additional LCDM simulation (LCDMa) 
uses a smaller force softening 77 = 10 /i~ 1 kpc and is evolved with 5000 time steps in order to 
check the possible effect of the force softening on the final dark matter distribution especially 
at small scales. As far as the shape of the virialized halos is concerned, we made sure 
that both simulations (LCDM and LCDMa) yield almost identical results. In what follows, 
therefore, we do not distinguish LCDM and LCDMa, and simply refer to them as LCDM. 
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2.2. Identification of dark halos in the cosmological simulations 

The Friends-Of- Friends (FOF) method is a widely used algorithm to identify dark mat- 
ter clumps in N-body data. The mean overdensity within the clumps is approximately 
proportional to b~ 3 , where b is the bonding length. It has been shown that the FOF clumps 
with b = 0.2d, where d = L/iV 1 / 3 is the mean separation of particles, approximately cor- 
respond to the virialized dark matter halos of the mean overdensity 180 (e.g., Davis et al. 
1985; Lacey & Cole 1994). On the other hand, a large fraction of the FOF clumps identi- 
fied with b = 0.2d are known to form a system of multiple virialized halos that are bridged 
via thin filaments (e.g., Suto, Cen & Ostriker 1992; Suginohara & Suto 1992; Jing & Fang 
1994, hereafter JF94). JF94 proposed to compute the overdensity around the local potential 
minima within each FOF clump to separate the virialized halos. While this can effectively 
achieve the goal, it is time-consuming to find the local potential minima (because there may 
be multiple minima within a single FOF clump). 

Here we propose to use an alternative method which works faster. The thin bridges 
connecting the halos identified with b = 0.2d can be effectively eliminated by reducing b. By 
trial and test, we found that the thin bridges almost disappear if we adopt b = O.ld. With 
this recipe, however, the resulting FOF clumps have a smaller size and a higher overdensity 
than those defined according to the spherical collapse model. Therefore our scheme should 
be interpreted to identify first the central parts or the substructures of the entire halo. Next, 
for each FOF clump of b = O.ld, we compute the gravitational potential of every member 
particle. The position of the particle of the minimum potential is defined as the center of 
the hosting halo. Then the spherical overdensity is computed around the halo center with 
increasing the radius, and the virial radius r vir is found when the overdensity reaches the 
value predicted in the spherical collapse model. Here we use the fitting formula of Bryan & 
Norman (1998) for the spatially-flat (fl(z) + X(z) = 1) models: 

A vir (^) = 3 f vir = 18tt 2 + 82 [Q(z) - 1] - 39 [Q(z) - if , (1) 

where p crit is the critical density of the Universe. Since our choice b = O.ld preferentially 
selects smaller clumps than those predicted in the spherical model, some fraction of such 
clumps turn out to be substructures within the virial radius of a larger halo defined in the 
above equation. If the virial spheres of more than one halos overlap, we simply retain the 
most massive clump and throw away the others from the final halo list. 
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2.3. High-resolution halo simulations 

Our cosmological simulations which we described above have a sufficient spatial resolu- 
tion to discuss the statistics concerning the halo shapes and the concentration of the density 
profile (§4 and 5) as was conducted by Jing (2000) in the framework of the spherical approx- 
imation. Actually except for a delicate problem of determining the slope of the central cusp 
at r <C 0.01r V i r , a larger simulation volume is more important than the higher resolution 
for the current purpose. Nevertheless we also use our higher-resolution halo simulations 
(Jing & Suto 2000; hereafter simply referred to halo simulations) to demonstrate that our 
triaxial modeling indeed provides a better description for halo profiles than the conventional 
spherical modeling (§3). 

These halos are simulated with about a million particles within their virial radii (see 
Table 1 of Jing & Suto 2000). For mass scales of clusters, groups, and galaxies, there are 
four halos, respectively, and thus twelve halos in total. They are simulated in the LCDM 
model except the fact that the fluctuation amplitude, <t 8 = 1 (Kitayama & Suto 1997), is 
a bit larger than our current choice as = 0.9. Another advantage of the halo simulations 
is that those halos are simulated with almost the equal number of particles independently 
of the mass of the halos, and thus the resolution relative to the virial radius and the halo 
mass is kept constant. This is not the case for the cosmological simulations in which massive 
halos would have a better resolution in terms of the number of particles involved. Thus the 
possible artificial effect due to the variable resolution is suppressed in the halo simulations. 

After Jing & Suto (2000) was published, we have completed runs of additional two 
halos with a galactic mass and with a group mass. Those new halos are referred to GX5 and 
GR5, respectively, according to our previous convention. While we add these two, we also 
eliminate two previous halos from the list of halos that we examine below; GR2 which shows 
a clear bi-modal structure, and GX1 which is seriously disrupted at z < 0.5 due to the tidal 
force of a nearby massive object. This is because the major purpose of analyzing the halo 
simulation catalogs is to check the validity of the triaxial modeling for typical halos. The 
fraction of those atypical halos is properly taken into account in the statistics drawn from 
the cosmological simulations. Thus the above replacement does not bias our conclusion. 



3. Modeling the non-spherical density profiles of dark matter halos 

In this section, we propose that a non-sphericity in the density profiles of dark halos 
is well described by a triaxial model on the basis of the detailed analysis of the halo sim- 
ulations. In fact, we demonstrate that the triaxial modeling significantly improves the fit 
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to the simulated profiles, at least for relatively relaxed halos, compared to the conventional 
spherical model. The statistical description including the probability distribution functions 
for axis ratios and the concentration parameters will be discussed in the next section using 
the cosmological simulations. 

3.1. Defining the iso-density surfaces inside individual halos 

The shapes of dark halos have been previously studied by many authors (e.g., Barnes & 
Efstathiou 1987; Warren et al. 1992; Jing, et al. 1995; Thomas et al. 1998), and it is already 
well known that they exhibit a significant amount of departure from spherical distribution. 
Those previous studies first compute the inertial tensor for each halo, and then compute 
the distribution of the axial ratios and the correlation of the direction of the principal axes. 
While this is a well-defined method to characterize the shape of halos in principle, we do not 
employ this for two reasons. 

First, this method assumes that we know in advance which particles belong to each 
halo. In reality this is not the case since we usually attempt to determine the member 
particles of a halo and its shape simultaneously. This is serious because the inertia tensor is 
sensitive to the outer boundary of the halo where the membership of particles is also difficult 
to define. Previous studies get around the problem by applying the procedure iteratively; 
first, all particles within a certain spherical radius from the center of halo are included to 
compute the inertial tensor and the resulting ellipsoidal configuration. Next, those particles 
outside the ellipsoid are thrown away from the member particles of the halo, and the inertia 
tensor is re-calculated. This procedure is repeated until the solution converges. While this 
method seems to work well in previous low-resolution N-body data, we were not able to 
obtain a good convergence in the case of our high-resolution halos. This is ascribed to the 
fact that our high-resolution halos retain a significant amount of substructures which have 
been artificially erased due to the overmerging effect in previous lower-resolution simulations. 
The iteration procedure is not stable in the presence of significant substructures especially 
at the boundary region of halos, since the inertial tensor is quite sensitive to them. 

Second, our main interest here is not simply to define the overall shape of halos, but 
to characterize the density profile. Therefore we would like to have a sequence of iso- 
density surfaces with different overdensities. The ellipsoidal surface obtained from the above 
procedure, even if it converges, is not related to those iso-density surfaces, and thus not so 
useful after all for our purpose here. 

With the above problems of the previous method in mind, we propose another approach 
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to find the iso-density surfaces. This begins with the computation of a local density at each 
particle's position. We adopt the smoothing kernel widely employed in the Smoothed Particle 
Hydrodynamics (SPH) method (e.g., Hernquist & Katz 1989): 



W(r,hi) = 



' ( 1 f 1 U>,< r< -lb,). (2) 



4 V K, 
otherwise 



where hi is the smoothing length for the i-th particle. We use 32 nearest neighbor particles 
to compute the local density pi, and hi is set to be a half of the radius of the sphere that 
contains those 32 neighbors. Using p iy we construct the iso-density surfaces corresponding 
to the 5 different thresholds: 

= A^Pcnt, (3) 

A (n) = 100 x 5"" 1 (n = 1 ~ 5). (4) 

In practice, we collect all particles satisfying 0.97pi™^ < pi < 1.03pi™' ) to define the n-th 
isodensity surface. The typical sizes (the mean radii) of those surfaces are 0.6, 0.4, 0.25 0.12 
and 0.06 times the virial radius of the halo, respectively. Note that pi^ is the local density, 
and thus the mean density of the halo inside the corresponding radius of p^ is generally 
much higher. 

Actually a straightforward application of equation (3) results in many small distinct 
regions with the identical density threshold inside an individual halo. This is again due to 
the presence of the strong substructures in the halo. Since we are interested in the isodensity 
surfaces which represent the overall density profile of the parent halo, we have to eliminate 
those small regions corresponding to the substructures. For this purpose, we again use the 
FOF technique but with a different bonding length from that we used when identifying the 
virialized halos. After some trial and error, we find that an adaptive (i.e., dependent on 
each isodensity value) bonding length of b n = 3(p^ /m p )~~ 1 ^ works well (c.f., Suto, Cen & 
Ostriker 1992). 



3.2. Triaxial model fits to the iso-density surfaces 

Figure 1 plots typical examples of the projected particle distributions within the iso- 
density surfaces for four different halos (CL3, GR1, GR5 and GX3) after particles in strong 
substructures are eliminated as described above. Those plots clearly suggest that the isoden- 
sity surfaces are typically approximated as triaxial ellipsoids. So we performed the following 
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triaxial fit to the iso-density surfaces with 5 different thresholds separately: 



R\p s ) 



X 2 Y 2 



Z 2 



(5) 



« 2 (Ps) b 2 ( Ps ) 
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The origin of the coordinates is always set at the center of mass of each surface, and the 
principal vectors a, b and c (a < b < c) are computed by diagonalizing the inertial tensor 
of particles in the surface (Fig. 2). The projected views of the corresponding fitted ellipsoids 
are shown at the bottom panels in Figure 1, which implies that the ellipsoid fitting is a good 
approximation (at least visually). 

Figure 3 plots the dependence of the axis ratios, a/c and b/c, on the isodensity threshold 
p s . Naturally each halo exhibits different behavior which may reflect the different merging 
history and/or tidal force field. Nevertheless, several systematic dependences are quite vis- 
ible. The halos of cluster mass generally have smaller axial ratios than those of galactic 
mass, implying that the halos of the galactic mass are rounder on average than those of 
cluster mass. This mass-dependence will be quantified with a large sample of halos from the 
cosmological simulations in §4. 

On the other hand, we also note that the axial ratios decrease with increasing the 
density; the iso-density surfaces become more elongated in the central region than in the 
outer region. The mean (with the error bar of the mean) of the axial ratios computed from 
the twelve halos are plotted in the right panels of Figure 3 (the symbols). The solid lines 
show the single power-law fit for the mean axis ratios: 



Figure 4 shows the degree of the alignment of the axis directions among isodensity 
surfaces at different densities (radii). We define On as the angle between the major axis of 
the isodensity surfaces and that of the = 2500 isodensity surface as shown in Figure 
2. Similarly, 2 2 is defined with respect to their middle axes. According to our definition, 

cos On = cos ^22 = 1 at p s /p C rit = -4 (3) = 2500. 

We find that the major axes align pretty well within a halo; for about 70% of the halos 
cos #11 at different radii is larger than 0.7. For about half of the sample, cos#n is larger than 
0.9. In a few cases (3 out of 12 halos), however, the alignment of the major axes is poor. 
When we check these halos individually (e.g. GR1), it turns out that b/c for the two halos is 
quite close to unity, indicating they are oblate halos with b ~ c and thus the direction of the 
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major axis is difficult to measure (if b = c, the direction of the major axis is arbitrary within 
a plane). Thus the apparent mis- alignment of their major axes is not meaningful. Only 
for the remaining one halo (GX3; Figure 1), the major axes of the outer and the innermost 
isodensity surfaces are indeed perpendicular to that at the middle. This is the real case that 
the major axes are significantly mis-aligned. 

The alignments of the middle axes show similar behavior: for most of the halos the 
degree of the alignment is satisfactory. For those which show significant misalignment of the 
middle axes, their a/b or b/c ratio is usually quite close to unity and the direction of the 
middle axes (and the minor or major axes) can be poorly determined at best. Only in the 
case like GX3, no simple ellipsoid description can be found, but this is fairly exceptional. 
The alignment seems slightly better for cluster-sized halos, but this would be simply because 
galactic halos are more spherical and thus the direction of the major axis is less accurate 
than that for cluster-sized halos. 



In the last subsection, we have seen that the isodensity ellipsoids at different radii are 
approximately aligned, and the axial ratios of the ellipsoids are nearly constant. These 
facts suggest the possibility that the internal density distribution within a halo can be 
approximated by a sequence of the concentric ellipsoids of a constant axis ratio. To show 
this to be an improved description over the conventional spherical description, we compute 
the quadrupole of the particle distribution within a spherical shell (Q s ) or an ellipsoid shell 
(Q e ). For a spherical shell, the positions of particles inside the shell can be described by 



y z = r cos 9 

with r being the (conventional) spherical radius. Similarly, the positions of the particles in 
an ellipsoidal shell can be described by 



where X, Y and Z axes are the principal vectors of the ellipsoidal shell, and a/c and b/c are 
the axis ratios. In the rest of the paper, we preferentially use the capital R to refer to the 
length of the major axis defined in the triaxial model. 



3.3. Triaxial versus spherical modeling of dark halos 




(8) 
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Then the quadrupole moments of the iso-density surfaces in the spherical and triaxial 
models, Q s and Q e , are computed as 

2 

- i, (io) 

2 

-1, (11) 

where the summation over j runs for all particles (N p ) in the iso-density surface, and Yi m 
is the spherical harmonics. If the spherical (triaxial) model is exact, Q s (Q e ) vanishes. 
Using these measures, we will show the extent to which the triaxial model indeed provides 
a significantly improved description for the simulated halos. 

In practice, we compute Q s {r) and Q e (R) for 5 shells of each halo at r — R — 0.65r v i r , 
0.35r vir , 0.2r vir ,0.12r vir , and 0.065r vir with the shell thickness Ar/r = AR/R = In 10 x 0.1 = 
0.23. Those shells are centered at the potential minimum of the halo. In the triaxial model, 
we assume that the shells have the same axis ratios and the same principal axis directions 
as measured from the isodensity surface at = 2500. Thus those shells do not necessarily 
correspond to the iso-density surfaces that we have discussed. Actually this treatment is 
important because otherwise the triaxial model (with more degrees of freedom) should always 
provide a better fit. Also this approximation is most likely what one would like to apply 
statistically to halos of visible objects, which would yield a practical and fair comparison 
between the spherical and triaxial models. 

In the top and the middle panels of Figure 5, we present the quadruple moments Q s 
and Q e for the twelve halos. The quadruple Q s for the spherical modeling increases nearly 
monotonically with the radius. The Q e for our ellipsoidal modeling (in its simplified version 
as described above) stays flat at R < 0.3r vir but increases with the radius at the larger radius. 
The ratio of the two quadruples is shown in the bottom panel of Figure 5 which indicates 
that our triaxial model, even simplified, fits the simulated halo profiles much better than the 
spherical model. For 10 out of 12 halos, the ratio, Q e (R)/Q s (r), is much smaller than 1 at all 
scales (r = R). Even for the remaining two halos (GX3 and GR5), the ratio exceeds unity a 
bit only at the largest radius, and the triaxial description shows a significant improvement 
over the spherical model. The ratio Q e {R)/Qs{i r ) seems to approach unity as r becomes 
closer to r vir . The reason might be that the subclustering is more prominent in the outskirt 
region than in the central region, since the subhalos are tidally stripped when they fall into 
the central region (see Fig. 1 of Jing & Suto (2000)). The strong subclustering makes both 
the spherical model and the ellipsoidal model difficult to accurately describe the complicated 
density distribution at r vir . But the figure also clearly shows that our triaxial model works 
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significantly better than the spherical model for r < 0.7r vir rs r 2 oo, i-e. almost the entire 
halo (the definition for r 2 oo will be given shortly). 



3.4. Density profiles in the triaxial model 



The next important task is to describe the density profiles in the triaxial model gener- 
alizing the previous results in the spherical approximation (NFW; Moore et al. 1998; Jing 
& Suto 2000; Klypin, et al. 2001). In the same spirit of the previous subsection, we do not 
perform the fit to the iso-density surfaces that we identified, but rather compute the mean 
density p{R) at the simplified triaxial shells (i.e., the same axis ratios and the axis directions 
for the entire halo as those measured from its isodensity surface at = 2500) within a 
thickness of AR/R — 0.12. 

Figure 6 plots the density profiles measured in this way for individual halos as a function 
of R. As in the spherical case, we adopt the following form: 



where R is a scale radius and 5 C is a characteristic density. Again following the definition 
of r 2 oo in the spherical model (within which the mean matter density is 200p cr i t ), we define 
a radius R e so that the mean matter density within the ellipsoid of the major axis radius R e 
is A e p crit with 



The non-trivial dependence of A e on the axis ratios in the above equation is chosen so that 
R e becomes a fixed fraction of the virial radius r vir (see Fig. 7 below). 

The best-fits to equation (12) for each halo are shown in Figure 6 for a — 1.5 (solid 
lines) and for a — 1.0 (dotted lines). Up to the resolution limit of the halo simulations 
(R/Re ~ 0.02), equation (12) yields a good fit both for a — 1 and for a — 1.5. If comparing 
the fits to the simulation data more carefully, however, a = 1 works better for the halos 
of cluster mass and a — 1.5 better for those of galactic mass, which is consistent with the 
finding of Jing & Suto (2000) in the spherical model (but see Fukushige & Makino 2001, 
for a different point of view). 



p{R) 5c 



(12) 



Pcrit (R/Ro) a (1 + R/Ro) 
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(13) 



We also introduce a concentration parameter in our triaxial model: 
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which is plotted in the upper panel of Figure 7 adopting a — 1.0 (crosses) and a — 1.5 (filled 
circles) in the fit. In what follows we will not address the issue related to the inner slope 
of the density profiles, and adopt a — 1. It should be noted, however, that our statistical 
results presented in the next section can be readily applied to the a = 1.5 case since the 
ratio c e (a = 1.5)/c e (a = 1) is always close to 1/2. 

Before moving to the statistical analysis of halos in the cosmological simulations, we 
note that the value of R e and thus that of c e are dependent on our specific definition of A e 
(eq.[13]). As the middle and bottom panels in Figure 7 indicate, both R e /r v[T and c e /c v ; r 
(where c v ; r is the ratio of the virial halo radius to the scale radius r s in the spherical model) 
remain constant 0.45) independently of the mass of the halos when we adopt equation 
(13) for A e . This property is quite useful in applying our results for a variety of theoretical 
predictions, since for a halo of given virial mass M vir : 

M vir = — r 3 ir A vir p crit , (15) 

the radius R e in our triaxial model is easily computed. It is also known that the 
function of the halo mass (NFW; Eke, Navarro & Steinmetz 2001) with the scatter described 
by the lognormal distribution function (Jing 2000; Bullock et al. 2001). Therefore, once the 
shape of a halo at a given mass is specified, the density profile of the halo is completely fixed. 
The statistical distribution function of the halo shape is discussed in the next section. 



4. Statistics of triaxial density profiles 

High-resolution halo simulations, like those used in the last section, is well suited for 
studying the detailed internal structures of individual halos, but the number of such halos is 
too small for a statistical description. Therefore we switch to the halo catalogs constructed 
from our cosmological simulations in order to study the probability distribution of the shape 
of halos. As emphasized in §2, the cosmological simulations employ N = 512 3 particles in 
a 100/i _1 Mpc box and thus the mass resolution is even better than that of individual halo 
simulations of the original NFW paper, for instance. 

We consider halos which contain more than 10 4 particles within the virial radius. The 
lower mass limits are 6.2 x 10 12 M Q and 2 x 10 13 M© in the LCDM and SCDM models, 
respectively. We also consider three epochs at redshifts z — 0, 0.5 and 1.0 to examine the 
time-dependence. At these redshifts, we have 2494, 2160, and 1534 halos in the LCDM 
model, and 1806, 879, and 263 halos in the SCDM model, respectively. 
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4.1. Probability distribution of axis ratios 

Following the prescription presented in the last section, we determine the halo shapes 
at the iso-density surfaces with = 2500. Since the typical radius of the surfaces is about 
0.3r vir , they are well resolved in our cosmological simulations; the force softening length is 
typically smaller by one order of magnitude. 

Left panels of Figures 8 and 9 present the ratio a/cof the minor axis a to the major axis 
c for halos from the cosmological simulations in the LCDM and SCDM models, respectively; 
solid, dotted and dashed histograms indicate the results for 10 4 < iV halo < 2 x 10 4 , 2 x 10 4 < 
-^haio < 6 x 10 4 , 6 x 10 4 < iV halo , where iV halo is the number of particles within the virial 
radius of each halo; in those Figures we use M 4 = iV ha i o /10 4 , and thus M 4 = 1 corresponds 
to M vir = 6.2 x lO^/i^M©, and 2.1 x lO^/i^M© for our LCDM and SCDM models). Top, 
middle, and bottom panels show the results at z — 0, 0.5, and 1.0. 

Two systematic trends are visible; the ratio is slightly larger for less massive halos, and 
decreases at higher redshifts. This motivates us to attempt the following empirical scaling 
for the axis ratio a/c: 

- (as (os/M^yWV-' 

r - = J Ur J ■ (16) 

where M* is the characteristic non-linear mass at z so that the rms top-hat smoothed over- 
density at the scale a(M*, z) is 5 C = 1.68. The M* at z = 0, 0.5, and 1.0 are 9.4 x 1O 12 /i _1 M , 
2.0 x 1O 12 /i _1 M , and 3.8 x 10 11 h~ 1 M e respectively for LCDM, and 8.5 x lQ l2 h^M Ql 
9.7 x 10 u h- 1 M Q , and 1.4 x HF/i^M© respectively for SCDM. 

Such scaled axis ratios f ac show a fairly universal distribution almost independently of 
the mass and the epoch (histograms in the the right panels of Figs. 8 and 9). The universal 
probability distribution function of the ratio f ac is well fitted to the following Gaussian: 

P(r ac ) dr» = -±- exp (- ^^ ) #-c (17) 

with <t s = 0.113. 

Next we decompose the joint probability distribution function of the axis ratios as 

p(a/c,b/c)d(a/c)d(b/c) = p(a/c)d(a/c) p(b/c\a/c)d(b/c) 

= p(a/c)d(a/c) p(a/b\a/c)d(a/b) (18) 

in terms of the conditional probability distribution functions, p(b/c\a/c) and p(a/b\a/c). 
The second equality holds because once a/c is fixed, the distribution of a/b is uniquely 
determined from that of b/c. Since we have shown that the distribution function p(a/c) 
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is well approximated by equations (16) and (17), we compute the conditional probability 
distribution p(a/b\a/c). Figures 10 and 11 plot the results for the LCDM and SCDM models, 
respectively. Different panels correspond to p(a/b\a/c) for different ranges of a/c. Solid, 
dotted, and dashed histograms indicate p(a/b\a/c) at z — 0, 0.5, and 1.0, respectively. 

The conditional functions appear to be insensitive to the redshift. In both cosmo logical 
models, they are accurately fitted to 



p(a/b\a/c) 



2(1 - r min ) 



1 - 



(2a/b 



1 ^mir 



(19) 



for a/b > r min , where r min = a/c for a/c > 0.5 and r min = 0.5 for a/c < 0.5. p(a/b\a/c) = 
for a/b < r min . 



4.2. Probability distribution of the concentration parameter 

We apply the triaxial density profile (eq.[12]) obtained in the halo simulations to the 
halo catalogs in the cosmological simulations. Considering the resolution limits, we adopt 
a = 1 and use the data points at rj < R e < r vir in the fit, where rj the force softening length 
(see §2). Since we do not address the innermost structures of the halos and rather focus on 
the value of the concentration parameter c e , this catalog has a sufficient resolution to yield 
an unbiased estimate (e.g., Jing 2000; Bullock et al. 2001; Eke, Navarro & Steinmetz 2001 
for discussion). As already found in the spherical model (Jing 2000), the distribution of c e 
in the triaxial model has a significant scatter even if the range of the halo mass is fairly 
specified reflecting the dependence of the merging history of the individual halo. 

The resulting probability distribution functions for c e are presented in Figure 12, which 
are well fitted by the lognormal distribution: 



p{c e )dc e = —=L — exp 

\/2 / K(T r 



(lnc e — lnc e ) 2 



d\nc e (20) 



with a dispersion of a Ce 0.3 both in the SCDM and LCDM models. The dispersion is 
slightly larger than the value estimated in the spherical model 0.2) for equilibrium halos, 
but is comparable to the value for all halos put together (Jing 2000). It should be noted here 
that despite the fact that the triaxial model is superior in describing the density distribution 
of halos, the scatter in c e is comparable to the scatter of concentrations in spherical profile 
fits, which probably means that the scatter originates from the halo merger histories rather 
than non-sphericity of the halos. 
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The probability distribution (eq.[20]) is completed by specifying the mean of the con- 
centration parameter c e . The result from our simulations is plotted in Figure 13 as a function 
of the halo mass at z — 0, 0.5 and 1.0. NFW proposed a semi-analytic fitting formula for the 
concentration c v ; r in the spherical model. 1 More recently Bullock et al. (2001) have shown 
that in their LCDM model (the parameters are similar to those of our LCDM model here) 
c v ; r of a given mass decreases with z as oc (1 + z)~ x . The redshift dependence is stronger 
than that predicted in the NFW recipe. Thus Bullock et al. (2001) have proposed another 
recipe which successfully describes the concentration c vir . Since we have already shown that 
the ratio, c e /c vir , is almost constant (Fig. 7), it is interesting to see if the formula of Bullock 
et al. (2001) also describes the behavior of c e in our triaxial model. 

In the LCDM model, we find that the redshift dependence of c e for a given mass is 
approximately oc (1 + z)~ l in good agreement with their result. In the SCDM model, 
however, our result of c e shows a stronger redshift dependence than their prediction. This 
is also likely that the fitting formula of Bullock et al. (2001) was designed to fit spherical 
concentrations, c v ; r , and is not applicable to the non-spherical case due to the evolution of 

Ce/ C v ; r . 

Following NFW and Bullock et al. (2001), we propose a new fitting formula for c e in 
the triaxial model: 

^^Vlfe) 1 ' (2i) 

In the above, z c is the collapse redshift of the halo of mass M (NFW): 

er fc W-Jc = 1 (22) 

y/2[a*(fM)-a*(M)] 2' 

where o~(M) is the rms top-hat mass variance at z = 0, 5 C = 1.68, 5 c (z) = 1.68/ D(z), D(z) 
is the linear growth factor, and / = 0.01. Solid lines in Figure 13 indicate the predictions 
of equation (21), implying that the formula describes our simulation results very accurately. 
In those plots, we adopt A e = 1.1 and 1.0 for the LCDM and SCDM models, respectively. 

We also made sure that the formula also agrees well with our halo simulations in the 
LCDM model, while the results appear 10 to 20 % higher (i.e., A e — 1.2 ~ 1.3) than those of 
the cosmological simulations (A e = 1.1). Considering both the typical 30% scatter in c e and 
the limited number of the high-resolution halos (12 in total), the above level of difference 



1 Originally NFW defined the concentration parameter as C200 = r 2 oo/r s , where r 2 oo is the radius within 
which the mean overdensity is 200p C rit- Their recipe, however, can be easily generalized to c v i r , since r20o/r v 
is almost constant for a given cosmology. 
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may not be interpreted so seriously at this point. In fact, the difference may be attributed 
partly to the fact that halos with significant substructures (like GR2) have been eliminated 
in the high-resolution halo samples (§2), while we have not attempted such a selection in 
the cosmological simulations. Indeed Jing (2000) has noted that halos in equilibrium are 
systematically more centrally concentrated than those with significant substructures. We 
also note that most previous studies including NFW and Eke, Navarro & Steinmetz (2001) 
have preferentially selected isolated halos in re-simulating with higher resolution, which 
would have less substructures and therefore have slightly higher concentration than average. 
If one is interested in halos in nearly equilibrium, the best-fit value of A e should become 1.3. 
Since c e /c vir remains constant (Figs. 7 and 15), the fitting formula (eq.[21]) can also be used 
for predicting c„ r in CDM models. 

Finally we have checked if the fitted values of R e and c e are dependent on the shapes 
of halos. Figure 14 presents the ratio of R e to the virial radius r vir as a function of the axis 
ratio a/b. Clearly R e /r v - u is independent of a/b and of the redshift (or equivalently the halo 
mass in unit of M+, see also Fig. 7), and approximately given by 0.45. Similarly, we find that 
Re/^vh is independent of b/c and a/c. On the other hand, the concentration parameter c e 
is slightly dependent on the halo shape. Figure 15 indicates that halos with smaller a/c are 
less centrally concentrated. 

In terms of the scaled axis ratio (a/c) sc (eq.[16]), the ratio of the mean concentration c e 
for a given f ac = (a/c) sc and the overall average c e (M, z) (eq.[21]) is well approximated by 



This fit is plotted in the solid line in Figure 15, which is in good agreement with the simulation 
data for different halo masses and both in the LCDM and SCDM models. 

In this section, we used the halos identified from cosmological simulations which have 
iV > 10 4 particles. The smallest halos are resolved much more poorly than the largest halos 
and the high-resolution halos (§3) that consist of ~ 10 6 particles. In order to make sure that 
our results are robust to the mass resolution, we repeat the same analysis of the axis ratios 
and the density profile by randomly selecting N = 10 particles from each of the twelve 
high-resolution halos of the previous section. The results are compared in Figure 16 with 
those obtained in the last section where we consider all particles. Both the axis ratios and 
the concentration of the randomly selected sample agree well with those of the original halo 
sample; the typical dispersion of the axis ratios between the two samples is ~ 10%, and the 
concentration c e of the randomly selected sample is slightly lower (~ 8%). This comparison 
indicates that the mass resolution does not affect our results in this section significantly. 




(23) 



-18- 



5. Summary and Discussion 

This paper has presented a triaxial modeling of the dark matter halo density profiles 
extensively on the basis of the combined analysis of the high-resolution halo simulations (12 
halos with N ~ 10 6 particles within their virial radius) and the large cosmological simulations 
(5 realizations with N = 512 3 particles in a 100/i _1 Mpc boxsize). In particular, we found 
that the universal density profile discovered by NFW in the spherical model can be also 
generalized to our triaxial model description. Our triaxial density profile is specified by the 
concentration parameter c e and the scaling radius Rq (or the virial radius R e in the triaxial 
modeling) as well as the axis ratios a/c and a/b. 

We have obtained several fitting formulae for those parameters which are of practical 
importance in exploring the theoretical and observational consequences of our triaxial model 
(in doing so we have adopted a = 1 since the precise value of the inner slope is difficult to 
reliably determine even with the resolution of the current simulations); 

• the mass and redshift dependence of the axis ratio, or equivalently the definition of the 
scaled axis ratio f ac = (a/c) sc : equation(16) 

• the probability distribution of the axis ratio p(r ac ) : equation(17) 

• the conditional probability distribution of the axis ratios p(a/b\a/c) : equation(19) 

• the mean value of the concentration parameter c e (M,z) : equation(21) 

• the dependence of the concentration parameter on the axis ratio f ac : equation(23) 

• the probability distribution of the concentration parameter p(c e ) : equation(20) 

Since c e /c V i r remains constant (Figs. 7 and 15), the fitting formula (eq. [21]) can also be used 
for predicting c V i r in CDM models. 

We have focused on the triaxial modeling and characterization of dark halos in the 
present paper, and plan to show specific applications elsewhere. Nevertheless it would be 
worthwhile to mention several important examples of the current model. 

The results of the paper are applicable in a fairly direct manner to the following three 
areas, (i) the weak and strong lens statistics. The comparison with the weak lensing ob- 
servations provides information of the degree of triaxiality of observed clusters, mainly at 
outer regions. In addition, the frequency of the lensing arc is known to be sensitive to the 
non-sphericity of the halo mass profile especially in the central regions (e.g., Bartelmann et 
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al. 1998; Meneghetti et al. 2000, 2001; Molikawa & Hattori 2001; Oguri 2002). (ii) predic- 
tions of the non-linear clustering of dark matter based on the halo model (e.g., Mo, Jing, 
& Bonier 1997; Ma & Fry 2000; Hamana et al. 2001; Kang et al. 2002). The high-order 
statistics of clustering, e.g. the three-point correlation and the bispectrum, should be quite 
sensitive to the non-sphericity, (iii) Dynamics of galactic satellites. Recently this has been 
argued to be very sensitive to the non-sphericity of the host halo (e.g., Ibata et al. 2001). 
The combination of those three approaches would even yield a direct test of the cold dark 
matter paradigm (Spergel & Steinhardt 2000; Yoshida et al. 2000). 

Of course the non-sphericity of dark matter halos is critical to understanding that of 
the gas density profile of clusters of galaxies. Since gas dynamics is characterized by the 
isotropic pressure tensor and does not directly follow the dark matter distribution in halos. 
In fact, most hydrodynamic simulations of galaxy clusters suggest that the gas distribution 
is generally rounder than that of dark matter. Nevertheless we would like to mention a 
couple of important examples where the non-sphericity in the gas distribution has crucial 
and observable consequences, (iv) the gas and temperature profiles of X-ray clusters. Al- 
most all previous analytical models for the X-ray profiles of galaxy clusters have adopted the 
spherical approximation perhaps due to the lack of any specific model for the non-sphericity. 
Since our triaxial model specifies the gravitational potential of the hosting halos, one may 
compute the gas or temperature profiles, with an additional assumption of the hydrostatic 
equilibrium for instance, as performed in the NFW model (e.g., Makino, Sasaki & Suto 1998; 
Suto, Sasaki & Makino 1998; Komatsu & Seljak 2001). If combined with the observed surface 
brightness distribution of clusters, one may in principle solve for the gas and temperature 
profiles simultaneously for a given non-spherical distribution of dark matter (Silk & White 
1978; Yoshikawa & Suto 1999; Zaroubi et al. 1998). (v) the systematic bias and statistical 
distribution of the Hubble constant estimated via the Sunyaev - Zel'dovich effect. In view 
of the on-going observational projects, it is of vital importance to re-evaluate the reliability 
of the estimates taking account of the non-sphericity effect of the clusters. With the above 
modeling of the gas and temperature profiles for individual clusters, one may discuss the sta- 
tistical properties of the estimates of the Hubble constant (e.g., Fox & Pen 2002), combining 
the extensive fitting formula for the probability distribution functions of the triaxial model 
parameters and the halo mass function (e.g., Sheth & Tormen 1999; Jenkins et al. 2001). 
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Table I. Model parameters for cosmological simulations with N = 512 3 in a 100/i 1 Mpc 

box. 



Model 


^0 


Ao 


0"8 


r 




7] [h 1 kpc] 


timesteps 


realizations 


LCDM 


0.3 


0.7 


0.9 


0.2 


6.2 x 10 8 


20 


1200 


2 


SCDM 


1.0 


0.0 


0.55 


0.5 


2.1 x 10 9 


20 


1200 


2 


LCDMa 


0.3 


0.7 


0.9 


0.2 


6.2 x 10 8 


10 


5000 


1 



Table 2. Properties of the new simulated halos in the LCDM model with f2 = 0.3, 



A = 0.7, h 


= 0.7, a 8 = 


1, and T = 


0.21 


identification number 


M[h~ l M & ] 


a N p b 


r vir [/i- 1 Mpc] c 


GX 5 


6.1 x 10 12 


945864 


0.373 


GR 5 


5.5 x 10 13 


644839 


0.776 



a Mass of the halo within its virial radius. 
b Number of particles within its virial radius. 
c the virial radius of the halo. 
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Fig. I. — Examples of projected particle distribution in four halos; a) CL3, b) GR1, c) GR5, 
and d) GX3. The size of each box is 2r vir of each halo. For each halo, particles in the 
isodensity shells with A = p s /p C rit = 100, 2500, and 6.25 x 10 4 are plotted on the xy, yz and 
zx planes (from left to right). The bottom panels show the triaxial fits to five isodensity 
surfaces projected on those planes. 
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Fig. 2. — A schematic illustration of the triaxial model for the isodensity surface of dark 
matter halos. The axis lengths are defined to be a < b < c, and 9 n {622) measures the angle 
between the longest (middle) axis of the iso-density surface with that of = 2500 (eq.[4]). 
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Fig. 3. — Axis ratios for the triaxial model fits to twelve halos. Left: results for individual 
halos. The dashed lines are for cluster halos, the dotted ones for group halos, and the solid 
lines for galactic halos. Right: symbols indicate the mean and its one-sigma error from the 
halo simulations, while the solid lines show the single power-law fit (eq.[6]). The upper and 
lower panels show a/c and b/c, respectively. 
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Fig. 4. — Degree of alignment of the directions of the ellipsoid axes. Left: results for 
individual halos. The dashed lines are for cluster halos, the dotted ones for group halos, and 
the solid lines for galactic halos. Right: symbols indicate the mean and its one-sigma error 
from the halo simulations. The upper and lower panels show for the major and middle axes 
respectively. 
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Fig. 5. — The quadrupole moments denned in the triaxial model (Q e ) and in the spherical 
model (Q s ) for five shells at radii from 0.05r vir to 0.65r vir . They are presented in the top 
two panels, and their ratio Q e /Q s is in the bottom panel. The dashed lines are for cluster 
halos, the dotted ones for group halos, and the solid lines for galactic halos. In the electronic 
edition, the cyan dotted line is for GR5 and the green solid one for GX3. 
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Fig. 6. — Radial density profiles in our triaxial model of the simulated halos of galaxy (left), 
group (middle), and cluster (right) masses. The solid and dotted curves represent fits to 
equation (12) with a = 1.5 and 1.0, respectively. For reference, we also show p(R) oc 
and R~ 15 in dashed and solid lines. The vertical dashed lines indicate the force softening 
length which corresponds to our resolution limit. For the illustrative purpose, the values of 
the halo densities are multiplied by 1, 10 _1 , 10 -2 , 10 -3 from top to bottom in each panel. 
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Fig. 7. — The fitting results of the triaxial model to twelve halos. a) the concentration 
parameter c e for a = 1 (crosses) and for a = 1.5 (filled circles); b) the ratio of c e to that of 
the spherical counterpart, c v i r , for a = 1; c) the ratio of R e to the virial radius r v - ir in the 
spherical model. 



-32 - 



=0.3, A =0.7, cr 8 =0.9 
raw data scaled 

6 I — i — i — i — | — i — i — i — | — i — i — i — | — i — i — i — | — i — i — i — I — i — i — i — | — i — i — i — | — i — i — i — | — i — i — i — | — i — i — r 




Fig. 8. — The distribution of the axis ratio a/c of the halos in the cosmological simulations 
of the LCDM model before (left) and after (right) the scaling described in the text. Top, 
middle and bottom panels correspond to z — 0, 0.5, and 1.0, respectively. Solid, dotted 
and dashed histograms indicate the results for halos that have the number of particles of 
M4 = (iVhaio/10 4 ) within the virial radius. The smooth solid curves in all the panels represent 
our fit (eq.[17]). 
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Fig. 9. — Same as Figure 8, except for the halos in the SCDM simulations. 
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Fig. 10. — The conditional distribution of the axis ratio a/b of the halos in the cosmological 
simulations of the LCDM model for a given range of a/c. Halos at different redshifts are 
represented with different lines as indicated in the bottom-right panel. The smooth solid 
curves in all the panels represent our fit (eq. [19]). 
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Fig. 12. — Distribution of the concentration c e of the halos in the LCDM (left) and in the 
SCDM (right) models for different halo mass M4 = (iVhaio/10 4 ). The smooth solid curves 
represent our log- normal fit (eq.[20]). 
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Fig. 13. — The mean of the concentration c e as a function of the virial mass in the LCDM 
and SCDM models. The solid curves represent our fitting formula (eq. [21]) at z — 0, 0.5, 
and 1.0 from top to bottom. The data point, labeled Moore99, is taken from the result of 
Moore et al. (1999), and is scaled according to our fitting formula (eq.[21]). 
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Fig. 14. — The ratio of R e to the virial radius for halos with different shapes in the LCDM 
(left) and SCDM (right) models at different redshifts. 



-39- 




Fig. 15. — The dependence of the ellipsoid concentration c e on the scaled axis ratio (a/c) sc . 
Different symbols denote the results of halos of different mass (-/V halo = 10 4 M 4 particles) 
in the LCDM and SCDM models. The smooth solid curve represents our fitting formula 
(eq.[23]). 
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Fig. 16. — The axis ratios a/c (crosses, upper panel) and b/c (triangles, upper panel), and 
the concentration c e (lower panel) of a subsample of halo particles randomly selected from 
the twelve high- resolution halos compared with those of the whole sample (abscissa). Each 
halo in the subsample has N p — 10 4 particles. 



